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Non-linear effects on driven oscillations are im- 
portant in many fields of physics, ranging from 
applied mechanics to optics. They are instrumen- 
tal for quantum applications [U [2]. A limitation 
is that the non-linearities known up to now are 
featureless functions of the number of photons N 
in the oscillator. Here we show that the non- 
linearities found in an oscillator where supercon- 
ducting inductance is subject to coherent phase- 
slips, are more interesting. They oscillate as a 
function of number of photons N with a period of 
the order of yN, which is the spread of the co- 
herent state. We prove that such non-linearities 
result in multiple metastable states encompassing 
few photons and study oscillatory dependence of 
various responses of the resonator. A phase-slip 
process in a superconducting wire is a topologi- 
cal fluctuation of the superconducting order pa- 
rameter whereby it reaches zero at certain time 
moment and in certain point of the wire [3]. 

Such a process results in a ±2ir change of the super- 
conducting phase difference between the ends of the wire; 
this produces a voltage pulse. Incoherent thermally- 
activated phase slips were shown to be responsible for 
residual resistance of the wire slightly below critical tem- 
perature |U [S]. At lower temperatures and in thinner 
wires phase slips are quantum fluctuations. Although 
resistance measurements indicate the quantum nature of 
the phase-slips |B] , they cannot prove a possible quantum 
coherence of phase-slip events. A set of other nanodc- 
vices [7J |8] have been proposed to verify the coherence 
experimentally. To facilitate this verification was the ini- 
tial motivation of our research. 

The inductance L of the wire brings about the induc- 
tive energy scale Er, — Qq/L 2 , where $o = irh/e is the 
flux quantum with Ti the Planck constant and e the elec- 
tron charge. It is usually assumed that experimental ob- 
servation of coherent quantum phase slips requires the 
phase slip amplitude E$ to be comparable with El [5]. 
The phase-slip amplitude Es depends exponentially on 
the wire parameters, so its value can hardly be predicted 
and it may be small. This is why it is important to be 
able to detect arbitrary small values of Es- Our idea is 
to use a driven oscillator. We prove that in this case the 
detectable values of Es are only limited by damping of 
the oscillator Es ~ TiT <C Tiuiq. There is an outburst 
of activity in applying super conducting oscillators for 
quantum manipulation purposes [5]. The inductance of 



such an oscillator may be either a thin superconducting 
wire [TUJ [TT] or a chain of Josephson junctions [TH IT5] . 
The multi-junction chains also exhibit phase slips and for 
our purposes are very similar to a wire. Typical experi- 
mental values for the main frequency and dissipation rate 
are uj q ~ 10 10 Hz and T ~ 10 5 Hz. 

This brings us to the system under consideration: the 
phase slip oscillator. The setup is shown in Fig. [TJi and 
the equivalent circuit in Fig. [TJa. For simplicity, we ne- 
glect the effects of the capacitance distribution along the 
wire attributing all the capacitance C to the "island" at 
the end of the wire. Without phase slips, the system is a 
linear LC one-mode oscillator. The a.c. component of the 




FIG. 1: a. Phase slip oscillator. A thin superconducting wire 
connects a lead and an island. The nearby gate electrode in- 
duces charge to the island. The wire is subject to coherent 
phase-slips, b. The inductance L of the wire and the capac- 
itance C of the island form an oscillator that can be excited 
with the gate voltage. The crossed diamond represents phase- 
slips of amplitude Es ■ c. The phase-slip induced correction to 
level n is proportional to the overlap integral of the oscillator 
wave functions (shown here for n = 5; 8). Oscillations of the 
wave functions give rise to oscillations inn. d. The phase-slip 
corrections to the levels of the oscillator for 7 = 0.92. Exact 
values (green dots) are fitted with Eq. [2] The photon distri- 
bution in several coherent states is plotted in red to illustrate 
the correspondence between the width of the coherent state 
and the period of oscillations. 
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FIG. 2: Phase slip induced correction 8\ of a driven oscillation 
versus detuning at F = 15I\ Real and imaginary parts are 
shown. Upper pane: \/~N versus detuning at the same force. 
It illustrates the oscillation period A(VjV) — tv/j. 



gate electrode excites the oscillator while the d.c. com- 
ponent induces constant charge q = CV g to the island. 
The oscillator is subject to small damping characterised 
by the energy loss rate r <C ujq, ujq — 1/VLC. The 
dynamics of the oscillator can be described by (f> - su- 
perconducting phase difference dropping along the wire. 
The phase <j) can take any value and is not restricted to 
the interval interval (—it, tt). Without phase slips the dy- 
namics are entirely linear with the inductive energy given 
by Ei J ((j)/27T) 2 . So to say, the phase does not know that 
is supposed to be quantized in units of 2ir. Consequently, 
the charge q on the island does not affect the dynamics. 
The phase slips shifting the phase by ±27r can be de- 
scribed by a Hamiltonian acting on the wave function of 
the system [8]: 



H s V(cp) = ^-^2 e ±i7Ui ' e ^((t) ± 2tt). 
2 ± 



(1) 



The effect of weak phase slips (E$ <C Tiojo) on the 
resonantly driven oscillator originates from the shifts 
E n = (n\Hs\n) to the otherwise equidistant levels of 
the oscillator. We immediately see from Eq. [T] that 
E n cx cos (irq/e), so the charge induced affects the quan- 
tum interference of phase slips with opposite shifts. Any 
effect of phase slips is thus periodic in gate voltage with 
a period 2e/C. The experimental observation of such 
dependence unambiguously identifies the quantum coher- 
ence of phase slips. As mentioned, we are more interested 
in the oscillatory dependence on the number of photons 
n. One envisages the origin of such a dependence from 
the fact that the energy shifts E n are proportional to 
overlaps of the oscillator wave functions shifted by ±2ir 
with respect to each other (see Fig. [lj.). The oscilla- 
tions of those wave functions which are in phase with a 
period Acf> (x are thus converted into oscillatory 

dependence on the photon number. 

The important parameter 7 = (2GqZ/tt)^ 1 / 2 mea- 



sures effective impedance of the oscillator Z = y/L/C, 
(with Gq = e 2 /irh) and defines the quantum fluctuations 
of phase cx 1 /^ 2 . Commonly, electrical resonators have 
7 ^> 1. However, superconducting wires provide signifi- 
cant kinetic inductance which may make 7 ~ 1 |121 113] . 
In this letter, we concentrate on experimentally accessible 
range 0.3 < 7 < 3. We have found that in this interval 
the shifts at any n can be sufficiently well approximated 
by the large n asymptotics: 



E„ 



2Es cos(irq/e)- 



s(2 7v ^-f) 

1/4 ' 



7T7 • n 



(2) 



This remarkable is the central point of our findings: 
it shows that the phase-slips add very unusual non- 
linearities to the oscillator. The period of oscillations 
reads An = 2-K^/nj^. At 7 ~ 1 it compares with the 
"width" (n) = %fn of a coherent state of the oscillator 
corresponding to the average number of bosons n. For 
larger 7, the phase-slip shift E n will make more oscil- 
lations at the scale of the coherent state width. This 
suppresses the effect of phase slips at large 7. 

The weak non-linearities induced by phase-slips are 
neither visible nor important unless the oscillator is res- 
onantly driven. Therefore we include the driving force 
2V ac (t) = V exp(i(u>o — u>)t) + h.c, with the detuning 
\lo\ <C ujq. It is convenient to normalise the driving force 
such that it enters the Hamiltonian as a combination 
%F(b + &t)/2, where b and fet arc the boson annihilation 
and creation operators. The force is then F — e^/irTiV . 

In the absence of phase-slips, this force F brings the 
oscillator into a coherent state with the amplitude given 
by [reference] 



X = (b) = 



F 



(3) 



A straightforward but involved perturbation theory (see 
Supplementary Material) gives the first order correction 
(oc Es) to this amplitude, valid at any 7 and fcsT, 

5X = 5(b) = -2£; s cos(7r <Z /e)£:(7) y^Ji(2 7 |A|), (4) 



with £(7) = exp(— ^-cotanh^). The factor J\ (the 
first order Bessel function) incorporates oscillations cor- 
responding to the oscillatory behaviour of the energy 
shifts. The exponential factor £(7) is best understood 
as the effect of averaging of these oscillations over the 
width of the coherent state. The first order correction is 
exponentially suppressed at high temperature and 7> 1. 
While this does not imply that all corrections vanish at 
7 > 1, we prefer to work at 7 ~ 1, where the exponent 
is ~ 1. We will also assume ksT <§; hu)Q. 

In the linear regime, F — > 0, the correction amounts 
to the frequency shift frui — > huj — 2Es^f 2 cos(irq/e)£(j). 
The correction becomes noticeable if it is of the order of 
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the line width, Eg — T, and can be revealed owing to the 
oscillatory dependence on gate voltage. 

However, the applicability of the linear regime is lim- 
ited to almost no photon excited, A -C 1. At larger 
driving, the correction slowly decays with increasing 
N = \X\ 2 . At N ^> 1, the correction becomes signifi- 
cant if Es > max(w, r)iV 3 / 4 . It is interesting to note 
that the oscillatory correction enhances the dependence 
on the detuning ui. This is why the correction becomes 
significant at much smaller Es, Es ~ max(w, T)^ 1 / 4 , if 
one concentrates on the derivative of the amplitude with 
respect to detuning, dX/du. We illustrate the scale of the 
correction and its oscillatory dependence on A in Fig. [4] 
and refer to Supplementary Material for details of these 
estimations. 

Let us go beyond perturbation theory, to the regime 
where the phase-slip correction becomes large, leading to 
qualitatively different physics. We present a comprehen- 
sive semiclassical analysis that captures the essence of 
the full quantum solution. 

In the semiclassical approximation, we replace N by 
a continuous variable. The non-linearities modify the 
detuning winEq.|w-4 W + dE(N)/MN, where E(N) 
is defined by Eq. (2j at 7 w 1. Squaring Eq. [3] yields a 
self-consistency equation for N at given F and ui [14] : 

F 2 

That suffices to make implicit plots N{F,u>). For com- 
mon non-linearities dE/dN is cx N. This gives either a 
single solution for N(uj) or three solutions corresponding 
to two metastable states. The oscillatory dependence on 
N changes this drastically. To elucidate, we plot N(u) 
in Fig. |] at fixed F = 1ST. At negligible E s , N(u) is 
a Lorentzian. Phase-slip corrections shift the curve hor- 
izontally, the magnitude of the shift oscillating with a 
period ~ yN. At sufficiently large Es this results in im- 
pressive characteristic "corkscrew" shape. At any given 
lo within the oscillator line width one finds a multitude of 
states that differ in N. About half of these states are sta- 
ble. We stress the tunability of this phase-slip oscillator: 
small changes of the driving force, detuning, or charge 
induced change the number of stable states, thereby en- 
abling easy manipulation of N. 

Generally, one expects fluctuation-induced switching 
between the available stable states. The semiclassical 
analysis does not account for that . Nor does it prove if a 
given metastable solution corresponds to a pure quantum 
state. It is also not clear if the semiclassical prediction 
for the metastable solutions works for the states with 
few photons. To understand this, we have performed 
numerical simulations with the full quantum equation 
Eq. [6] for density matrix. For illustration, we set u> = 
and Es to a moderate value of 6 TiT. We initialise the 
density matrix to vacuum, |0)(0|, and compute its time- 
dependence while making a linear sweep of F from to 
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FIG. 3: Multiple stability in phase-slip oscillator. We plot 
the number of photons JV versus detuning at F — 15T as 
predicted by semiclassical approximation (Eq. [5|. In the ab- 
sence of phase-slips, N(ll>) is a Lorentzian (green line). While 
at Es = 5 hF the deviation from Lorentzian is barely visi- 
ble, the enhanced dependence on dN/dui brings the oscillator 
close to the threshold of bistability, as it is shown in inset. 
At Es = 20 HT there are already multiple intervals of u) with 
two stable configurations. At Es = 200 TiV the curve takes a 
pronounced corkscrew shape. There is a multitude of stable 
configurations in all the range of the plot. 

5.5 r and back. From semiclassics we expect up to 3 
metastable states in this force interval. Plotting (n) ver- 
sus F for different sweep durations T gives a series of 
curves with evident hysteresis (see Fig. |4^) . Generally, 
one expects the relaxation time of the density matrix to 
be of the order of 1/r. Remarkably, a noticeable hys- 
teresis persists even at time intervals 4 • 10 4 r _1 . This 
clearly indicates an exponentially long life-time of the 
metastable states even for few photons. We hypothe- 
sise that the oscillator spends most of the time in one 
of such states while rare switching between these states 
result in equilibration of the probabilities to be in these 
states. Such equilibration occurs at the time scale cor- 
responding to the slowest switching rate. To prove this 
illustratively, we have computed the equilibrium density 
matrix at F = 4.85 V and expanded it into a sum over 
orthogonal quantum states. We have found that the den- 
sity matrix is mainly contributed by three pure states: 
one "dark" state ~ |0) and two coherent-like state cen- 
tred around 5.5 and 16.3 photons respectively with prob- 
abilities 0.46,0.25 and 0.15. The remaining probability 
corresponds to "excited" states that have nodes at posi- 
tions of the coherent-like states centred at 5.6, 16.3. The 
relaxation time that characterises the slow switching is 
300 r^ 1 at this value of F. About 4000 photons are 
absorbed and emitted in the oscillator during this time 
interval; this proves the extraordinary robustness of the 
states involved. 

To conclude, we have investigated the effect of non- 
linearities produced in a superconducting resonator by 
coherent phase-slips. These non-linearities are very dis- 
tinct from those previously known owing to the oscil- 
latory dependence on number of photons with a period 
~ y/~N. We have demonstrated that at semiclassical level 
the non-linearities result in a multitude of metastable 
states. The position and number of these metastable 
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FIG. 4: a. Hysteresis in phase-slip oscillator. Crossed curve: 
semiclassical prediction for N(F) at ui = 0. The solid curves 
give the result of slow sweep of F from to 5.5F. The sweep 



duration takes values 10 , 2 



10 3 , 10 4 and 4 



10 4 F" 1 from 



lowest to uppermost curve. The pronounced hysteresis indi- 
cates exponentially long life-times of the metastable states. 
Black circles indicate the pure states that contribute to the 
equilibrium density matrix at F = 4.85F (in detail in Fig.b). 
b. Pure states. We show diagonal elements of the equilib- 
rium density matrix at F = 4.85r. The pulses of different 
colour show contributions to diagonal elements from three 
pure states of highest probability. We see a "dark state" with 
predominantly zero number of photons and two coherent-like 
states centred around 5.5 and 16.3 photons. 



Here E(n) = (n\H\n) are the energy shifts induced by 
the phase-slips. These shifts are expressed through the 
hypergeometric function iF\. 

E n = 2E S cos(7r g /e) • exp (-772)^1 [-n, 1, -y 2 ]. 



Large-n asymptotics are given by Eq. [2 
slips (E n = 0), the solution of Eq. 
intervals reads 



Without phase- 
(16} for long time 



p(t — > 00) = exp 



(tfb-\tf + X*b + \X\ 2 



where A = (b) is defined by Eq. [3] To compute the first 
order correction SX, we develop a rather involved pertur- 
bation theory for the density matrix that is outlined in 
the Supplementary Material, which leads to Eq. |4j 

For our numerical calculations, we solve the evolution 
equation Eq. [6] by 4th order Runge-Kutta method with 
a fine time step(~ 0.02 1/T). Typically, we take into 
account N = 50 states. Since we investigate slow dy- 
namics induced by rare switching, the simulation runs 
take rather long time (w lOh for T w 4 ■ 10 3 1/T). 



states can easily be tuned by changing the driving force. 
At quantum level, we have demonstrated that there is a 
single quantum state corresponding to the semiclassical 
metastable states. These states are robust with an expo- 
nentially long switching time although they encompass 
only few photons. These features of the phase-slip os- 
cillator may be applied for a wide range of applications, 
as ultra-sensitive measurements, quantum manipulation 
and an unambiguous experimental verification of coher- 
ent phase-slips. 

We are grateful to J.E. Mooij, members of his team and 
to A. Ustinov for useful discussions. This work is sup- 
ported by the 'Stichting voor Fundamenteel Onderzoek 
der Materie (FOM)' and the 'Nederlandse Organisatie 
voor Wetenschappelijk Onderzoek (NWO)'. 

Methods 

We start with the Hamiltonian of the driven oscillator, 

H = huj b J< b + ReiFe^ ^}^^ + H s , 

where the phase-slip term Hg is given by Eq. [T] 

We implement the rotating-wave approximation to ar- 
rive to the equation for density matrix p : 

^- = - l -[H R ,p] + T (bptf -l( b i b p + ptfb)), (6) 

where the terms including the dissipation T in Eq. [6] are 
of conventional form [14] (assuming /c^T <C Tiloq) and 

H R = E(tfb) + Fbf+ 2 F * b +^b. 
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PERTURBATION THEORY FOR DENSITY MATRIX 



Here we present a more detailed account of how we have calculated the first order correction to the oscillator 
amplitude (Eq. |4| using perturbation theory for the density matrix that obeys Eq. [6j We rewrite it in compact form 

using a super-operator E that accounts for dissipation, detuning and the driving force: 
The stationary solution in the absence of phase-slips (Es = 0) reads 

hu>o 



exp 



k B T 



(tfb- A&t + A*6+|A| 2 ) 



where A = (b) is defined by Eq. [3] The phase-slip Hamiltonian E(tfb), Eq. [T] is treated as a perturbation. The first 
order correction to the density matrix satisfies the following condition 

d ^L = ± p W-i[E(tfb), Poo }. (7) 



This equation can be formally solved by making use of the superoperator resolvent R = (d/dt — E) 



-1 . 



P {1) {t) = - 1 - f dt'h{t-t')[E{tfb), Poo ]. 

n J -co 

To avoid explicit evaluation of this cumbersome resolvent we re-write it as an integral over a time-dependent operator 

poo ■ A 

P (1) (t) = J o drp(r), p(r) = --R(r)[E(bh), Poo ]. 

The thus defined p(r) obeys the time evolution equation 

dp(r) i 



dr h 

with the initial condition 



P,p(r)], (8) 



p{r = 0) = -^[E{b% Poo ]. 

Let us now represent the quantity E{b^b) via an integral that singles out the diagonal elements of the original phase-slip 
operator: 



Jo 



27T J 

C W e ixb f b Hse -i X b i b 

2tt 
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This operator consists of two shift operators that are expressed in terms of boson operators as e ±l i( bi + b ) . Therefore 



E(tfb) = E s 



2tt 



2tt 



Let us concentrate on the first order correction to the amplitude due to the 2tt phase-slip (corresponding to terms 
oc e l7Tq / e ). It is expressed as 



5(b) =Tr(6p«) = j dTTt[bp(r)] = -j.E s e i ^/ e j dr j d X Tr(b(p + (r, X ) 



P-(r, X ))). 



Here, /5±(r,x) satisfies Eq.^with initial conditions p+(0,x) = e^ (b+e,x+be ' x )p 0o and p-(0,x) = p^e^ elx+he * x \ 
We concentrate on the quantities £± defined as 

t±(r, X ) = ^ T f , T( X , 7 )=Tr(p ± (0 )X )). 
T(x,7) 

They obey the equation <9£±/<9t = — iw£± — iF — £±T/2, the same equation as the time dependent amplitude of the 
oscillations obeys. This yields 



with 



r(7) 



£fc(r,x) = r (7) + (Ci(x;7) - r)e- ( 5 +i " )T , 

Tr(6p±(0, X )) 



and £±(0,x) 



r(x,7) 



The integral over the time variable r gives a factor of ^ + iuj in the denominator. Let us also note that 

_ _ Tr([b, e' 7(bte ' x+be " x) ] ■ p^) 
^+ " T(x, 7 ) 7 ■ 

Combining this with the contribution of the — 2n phase-slips, we obtain 



6(b) 



~ + ILU Jo 2tt 



3 < X [ e We T ( X)7 )_ e -W«T( x ,-7)]. 



One computes the traces involved by making use of the expression for p, 

t 
2 



3^ q/e T( X , 7) - e-" q/e T( X , -7) = cos(7r g /e) exp 



cth (|) + i 7 Re (e* x A*) 



where a = ji^- To complete the calculation, we shift x with the phase ( of A* to obtain 



,,v E s hcos{-Kq e) 

(b) = ^— exp 

2+ lUJ 

Integration over x yields the final result, 



-^cth /Q 



^7 / ^2cos(x)exp[2i7|A|cos(x)]. 



2iE s cos(TTq/e)/h ( 7 2 /a\\ |A| 
<&> = r-r- exp (- y cth (-j j F7 Ji(27|A|). 

This is the oscillating function of g and A Eq. [4] discussed in the main text. 



ESTIMATIONS OF THE RELATIVE MAGNITUDE OF THE CORRECTION 

In this section we show how we have obtained the estimates for the phase-slip amplitude that causes the significant 
correction to the amplitude of driven oscillations. We use Eq.[4]and consider several limits. In the first limit, 7|A| <C 1, 
the first order Bessel function is expanded as J\{x) w f . Then the correction becomes significant if 



E S 



> max(w, r/2) 
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In the opposite limit, -y | A| 3> 1, the Bessel function is approximated by the cosine: Jiix) ~ y — cos(a;), and therefore 
the correction becomes significant if 

> max(a;,r/2) 3/4 

CJS ~ 7 l/2 

where N is the average number of excited photons. The correction relatively decreases with increasing N and in the 
limit N ^> 1 it oscillates. 

As mentioned in the main text, the oscillations enhance the frequency dependence of the correction. This feature 
might be used for experimental detection of the correction. To prove this, let us consider the derivative of the 
correction with respect to u> and compare it with the same derivative of the zero-order amplitude: 

dSX ,d\X\ E s r-— ,/n . _ . .. 

By analysing this we understand that the correction becomes significant if 

max(o;,r/2) 1/4 

E s <> N 

We see that this improves the estimation of the significant phase-slip amplitude by a factor of \/N. 

To summarise, we have three estimations for the significant phase-slip amplitudes. In increasing order, at Es — T, 
the correction is significant in the linear response regime iV -C 1. At Es ~ TN 1 / 4 one sees the significant effect on the 
frequency dependence of the oscillation amplitude (b). At Es ~ IW 3 / 4 , the overall change of the magnitude becomes 
dominant. 

We illustrate these conclusions with plots. We plot the relative magnitude of the correction versus detuning in all 
the figures, at a fixed value of the a.c. driving force F = 15T. This corresponds to N = 900 excited photons at 
zero detuning. At larger detunings this number decreases reaching TV ~ 1 at cu ~ F. Fig. [5] represents the real and 
imaginary part of the correction 6(b) /(b) from Eq. [4] The real part of the correction reaches its maximum at U) « F, 
as seen in Fig. [5j while the number of photons is small N ~ 1. Therefore one does not need to excite the oscillator 
in higher states to get large enough signal. At small frequencies it oscillates strongly with smaller amplitude, while 
at larger frequencies it decays to zero as l/u>. The imaginary part oscillates also strongly at small frequencies, but 
it goes to zero much faster, 1/w 2 . We also notice that the real part is odd in frequency, while the imaginary part is 
even. Their oscillations are in phase for uj > and opposite in phase at uj < 0. Fig. [6] shows the plot of the ratio 
of derivatives with respect to the frequency. We see how this enhances the relative magnitude of the corrections at 
small detunings (large TV). 
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